Correlation energies by the generator coordinate method: 
computational aspects for quadrupolar deformations 
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We investigate truncation schemes to reduce the computational cost of calculating correlations by 
the generator coordinate method based on mean-field wave functions. As our test nuclei, we take 
examples for which accurate calculations are available. These include a strongly deformed nucleus, 
156 Sm, a nucleus with strong pairing, 120 Sn, the krypton isotope chain which contains examples of 
soft deformations, and the lead isotope chain which includes the doubly magic 208 Pb. We find that 
the Gaussian overlap approximation for angular momentum projection is effective and reduces the 
computational cost by an order of magnitude. Cost savings in the deformation degrees of freedom 
are harder to realize. A straightforward Gaussian overlap approximation can be applied rather 
reliably to angular-momentum projected states based on configuration sets having the same sign 
deformation (prolate or oblate), but matrix elements between prolate and oblate deformations must 
be treated with more care. We propose a two-dimensional GOA using a triangulation procedure to 
treat the general case with both kinds of deformation. With the computational gains from these 
approximations, it should be feasible to carry out a systematic calculation of correlation energies 
for the nuclear mass table. 

PACS numbers: 21.60.Jz, 21.10.Dr 
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I. INTRODUCTION 

Much progress has been made in developing a fully 
microscopic method to determine nuclear binding ener- 
gies. The most successful approach in this framework [j 
starts from an energy functional of the nuclear orbitals 
0, which is minimized by solving mean- field equations 
for the orbitals. However, to achieve a level of accuracy 
around 700 keV on the more than 2000 known masses, 
some correlation energies beyond what can be subsumed 
within the mean-field functionals have to be introduced. 
Unfortunately, the well-established microscopic methods 
which could be applied are far too costly to be used for 
such a large number of nuclei. As a consequence, in 
the systematics of binding energies, the effects of cor- 
relations are estimated with phenomenological methods 
whose range of validity is not apparent. Our goal is find 
microscopically founded yet easy-to-implement methods 
allowing an accurate calculation of these correlations. 

There are several systematic approaches to correlation 
energies. Among them, the generator coordinate method 
is especially promising, mainly because of two attrac- 
tive features: it is a fully variational method and it is not 
limited to small-amplitude collective motion. It has been 
used successfully by the Paris-Brussels collaboration for 
calculating spectroscopic properties of nuclei (for repre- 
sentative applications, see Refs. SSIlllHUSGil), 

and we believe it has promise for a global microscopic 
theory to systematically calculate nuclear masses. 

The most important fluctuations beyond the mean field 
are those associated with pairing fluctuations and with 
quadrupolar deformations. We only discuss the latter in 
this work; the pairing fluctuations are treated by num- 
ber projection as in the above cited GCM studies. At 



the mean-field level, quadrupolar deformations are in- 
troduced by using a quadrupolar constraining field to 
generate deformed configurations. To go beyond the 
mean field, one mixes configurations to obtain an addi- 
tional correlation energy. This can happen in two ways. 
First, whenever the mean-field minimum is deformed, it 
can be interpreted as the intrinsic state of a rotational 
band. The lowest mean-field energy then corresponds 
to a weighted average over the members of the band, 
and one gains correlation energy by projecting onto the 
ground state of the band. Second, fluctuations of the 
quadrupole deformation about the value of the mean- 
field energy minimum also contribute to the correlation 
energy. In fact, when that minimum is spherical, they 
are the only one present. Both kind of correlations can 
be introduced simultaneously by projecting good angular 
momentum from the intrinsic states and by mixing con- 
figurations around the mean-field energy minimum with 
the generator coordinate method with the quadrupolar 
deformation as generator coordinate. Note that because 
the resulting wave function has amplitudes for deformed 
as well as spherical configurations, the distinction be- 
tween these two correlations becomes blurred. 

In previous works 0, IU , the procedure to in- 
troduce these correlations consists of three steps. First 
one constructs mean-field wave functions for a set of con- 
figurations that are defined by a constraining quadrupo- 
lar field. These configurations may be labeled by the 
Bohr-Mottelson variables of two deformation parameters 
(3 and 7. From them, additional states are generated by 
rotations with Euler angles SI = (<fr, 8, x) from which are 
projected states of good angular momentum. The pro- 
jected states corresponding to different values of (3 and 7 
are then mixed by the GCM, which means numerically 
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FIG. 1: Overview of potential energy landscapes for selected nuclei. See text for details. 



solving the Hill-Whcclcr equations in the discrete basis 
of the projected configuration set. The computational 
issue in such calculations is in the selection of the Euler 
angles and the deformation parameters that are neces- 
sary to achieve a sufficient accuracy. Up to now, a fairly 
large number of active angles (~ 12) were used for each 
deformation. We will show that this can be drastically 
simplified by using a suitable Gaussian overlap approxi- 
mation, particularly the topological Gaussian overlap ap- 
proximation advocated in ref . |ll| . 

Figure Ogives an overview of the energetics associated 
with the (3 degree of freedom in selected nuclei. We have 
singled out 208 Pb as an example of a doubly magic nu- 
cleus, 86 Kr as a soft spherical vibrator, 186 Pb as a soft 
nucleus with several coexisting near-degenerate minima, 
and 156 Sm as a well-deformed nucleus. The upper panel 
shows the energy surface, measuring the energy with re- 
spect to that of the GCM ground state. The dashed 
curve shows the energies of the states obtained by N 
and Z number projection from the HF+BCS configura- 
tions. The configurations are labeled by their intrinsic 
quadrupole moments 1 
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with ro = 1.2 A 1 / 3 fm. In the figure, the states are con- 
nected by the line. One sees that the ground state is 
spherical except for 156 Sm, which has a strong prolate 
deformation. The solid curves show the results of angu- 
lar momentum projection onto J = 0, labeling the states 



1 The negative values of fj occur for oblate axially symmetric de- 
formations; in the Bohr-Mottelson parameterization they corre- 
spond to /3 = |/3 2 | and 7 = 60°, 180°, or 300°, depending on the 
choice of symmetry axis. 



by the deformation of the state from which they were 
projected. One sees that there is a net gain of energy 
using a deformed intrinsic configuration for all the nu- 
clei including 208 Pb. The final energy, obtained by the 
Hill- Wheeler equations using the Z, N, J-projected con- 
figurations, is at zero energy on the graph. The dot in- 
dicates the average deformation of the mean-field states 
within the Hill- Wheeler ground state. The lower panel 
shows the relative amplitudes of the different configu- 
rations in the ground-state wave function. For 156 Sm, 
the function has a single peak and can easily be inter- 
preted as the zero-point motion of the p- vibration. For 
the other nuclei, there is a strong mixing of oblate and 
prolate configurations. As we will see below, this makes a 
simplified treatment for spherical nuclei somewhat more 
complicated than for the fully deformed nuclei. 

The GOA has often been used in the past to derive the 
ingredients of a collective Bohr Hamiltonian in the spirit 
of the work of Kumar and Baranger (12T. It is supple- 
mented by several other approximations [2, Il3| . The re- 
sulting Hamiltonian is five-dimensional, three dimensions 
corresponding to rotation and two to quadrupolar vibra- 
tions. This approach has been followed in ref. [14] and 
[15] to improve on the mean field theory of the Skyrme 
and Gogny interactions, respectively. The Lublin group 
has derived an alternative scheme starting from a Nils- 
son Hamiltonian 0, 0, 0] , which allows one to intro- 
duce several collective degrees of freedom. In the present 
work, the GOA is a numerical approximation, in prin- 
ciple under control, used to evaluate matrix elements of 
the norm and Hamiltonian kernels of the generator co- 
ordinate method, including projection on good angular 
momentum. 

Most previous calculations with the angular- 
momentum projected GCM have assumed axial 
symmetry, i.e. taking only configurations with 7 a 
multiple of 7r/6, and we will make the same assumption 
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here. In practive, the optimal spacing of configurations 
as a function of [3 remains to be determined, and also 
which off-diagonal matrix elements need to be calculated 
explicitly rather than estimated by the GOA. The first 
question was studied in ref. |l9j for unprojected matrix 
elements of the nucleus 194 Hg. We will see that the same 
considerations apply to the projected matrix elements, 
provided that the special characteristics of the [3 = 
singular point are taken into account. 

In the next section, we examine several methods for the 
approximate angular projection. Our target is to design 
a procedure leading to an accuracy of a few hundred keV, 
the level required to improve the best present mean-field 
based approach to binding energies. We find that a two- 
point approximation gives acceptable accuracy in all but 
rare circumstances. In Section IV we discuss the mesh of 
GCM configurations entering the Hill- Wheeler equation, 
and we propose a prescription for simplifying the space 
without losing the desired accuracy. 

Let us now briefly describe our procedure for dealing 
with pairing fluctuations. The mean-field calculations 
are performed in the Hartree-Fock+BCS approximation 
using the Skyrme interaction for the particle-hole channel 
field and a separate zero-range interaction for the pair- 
ing channel. The mean-field interaction is taken to be 
the SLy6 parameterization in the present work, the same 
as was used in ref. The pairing interaction is taken 
as a density-dependent delta function as in ref. |2Cj : it 
is used with 5 MeV cut-off of the orbital space above 
and below the Fermi energy. We use the same strength 
of the pairing as in ref. [2(|, (-1250 MeV fm 3 ), except 
for the Kr isotopes, where it taken as -1000 MeV fm 3 . 
To avoid the phase transition of the BCS theory at nu- 
clei with low densities of state at the Fermi surface, the 
Lipkin-Nogami prescription is used to generate the BCS 
wave functions, which are then projected onto fixed num- 
bers of protons and neutrons. The last step is essential to 
have nuclide-specific predictions. We always use number- 
projected wave functions in this work, and will refer to 
them as the intrinsic configurations of the GCM to dis- 
tinguish them from the BCS configurations that do not 
have definite particle numbers. 



II. ANGULAR PROJECTION 

The angular projection of a deformed wave function is 
carried out by calculating the integral 

I (q,q')=2 f dcos(6)(q\R e \q') (2) 
Jo 

and similar expressions for the matrix elements of the 
Hamiltonian and other operators. Here we have used the 
axial symmetry and the symmetry with respect to 180° 
rotations, (q\Rg\q') = to reduce the integra- 

tion over Euler angles to the above single integral. These 
properties can be used if the left and right wave functions 



have the same symmetry axes. We shall define our in- 
trinsic configurations accordingly, taking the symmetry 
axis along z. This requires that the oblate wave functions 
correspond to 7 = 180° instead of the usual 60°. 

In ref. |9( the integration was carried out on a Gauss- 
Legendre mesh with enough points to thoroughly sample 
the integrand in the range that it is nonzero. This is 
done in such a way that using the estimate for the over- 
lap given in Ref. [21j , the total number of abscissas in the 
interval [—1,1] is chosen in such a way that at least 12 
abscissas can be expected to give overlaps that are larger 
than 10~ 8 times the overlap at 9 = 0. The calculation 
of the integrand is started at 6 = and stopped when 
the calculated overlap has indeed fallen below this value. 
In this way, the choice of the Gauss-Legendre mesh is 
self-adjusting to the structure of the integrand, which 
scales with the dispersion of the angular momentum and 
is therefore nearly constant for states with small defor- 
mation and sharply peaked at 9 = for well-deformed 
ones. 

As discussed in the Appendix, the rotation operation is 
a rather costly computational task. The Gaussian over- 
lap approximation reduces the number of points to one 
(besides the identity) , making it attractive for global ap- 
plications. There are a number of ways that the approxi- 
mation can be applied. We shall write a generalized GOA 
for the overlap in the form 

(q\ReW) = (q\q')e- c ^ F W (3) 

where F has a fixed functional dependence on 9. Since 
there is only parameter in the expression, namely c(q, q'), 
only one rotation and many-particle overlap has to be 
carried out to determine it. The integral over the Gaus- 
sian is then a trivial computational task. 

The naive GOA takes F{9) = &\ which was proposed 
in early papers on the subject [22t |23| . Another form 
suggested in ref. [24[ is F = 1 — cos(0). We will call 
this the improved GOA. It is designed for cranking wave 
functions which break the time-reversal invariance. In 
our case here, the wave functions are invariant under 
time reversal and have additional rotational symmetries 
as well. With quadrupolar amplitudes to be rotated, a 
form having the correct limit for small deformations is 

m 

(q\Re\q') = exp { - ^c m , m / [6 mim i - P^ m ,(Q)] } 

where T> 2 is Wigner's rotation matrix for J = 2. Hagino 
et al. [ll| have considered this form for the case of axial 
symmetry, taking F as 

F = sin 2 -Vl . (4) 

Testing the approximation in a three-level Lipkin model, 
they found an important gain in accuracy using the 
sin 2 (9) dependence, which they call "topological GOA" 
(topGOA). They also tested the approximation on the 
Interacting Boson Model and found similar improvement 

n 
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FIG. 2: The overlap (q\R 9 \q) for the nucleus 208 Pb at q = 500 
fm 2 . The overlaps of the particle-number projected BCS in- 
trinsic states have been renormalized to one at 6 — 0. The cir- 
cles show the full calculation with the rotated many-particle 
wave functions. The fit by the topGOA is shown by the solid 
line and the fit by the improved GOA is shown as the dashed 
line. The 8 = point and the point shown as the filled circle 
were used to make the fit. 



We show in Fig. |21 the overlap (q\Rg\q) as a function 
of angle for a case of moderate deformation, 208 Pb at a 
deformation of -5 b (02 = —0.063). The points were cal- 
culated with the full many-particle wave functions. The 
GOA is constructed using the overlap at cos(0) = 0.75, 
where it is close to 0.5. The curves show the fits with 
the improved GOA and with the topGOA. While the 
improved GOA loses accuracy away from the peak, the 
topGOA gives a reasonable fit throughout. The projected 
integral Iq is off by 18% for the improved GOA and only 
1.3% for the topGOA. 

Next we consider the Hamiltonian matrix elements and 
their projection. Part of the GOA scheme as is usually 
applied assumes that the Hamiltonian matrix elements 
can be written as a product of the overlap times a factor 
that is quadratic in the difference of generator coordi- 
nates. For all versions of the GOA this may be expressed 



(q\HRg\q' 



-c(q,q')F(8) 



ho(q, q') 



h 2 (q,q')F(e) 
(5) 

The values of ho and hi may be determined from the two 
Hamiltonian matrix elements at 6 = and at the value 
used before to determine c. In Fig. [21 we show how well 
this works in the 208 Pb example. Plotted as circles are 
the ratios (q\HR$\q) /(q\R$\q) calculated from the many- 
particle wave functions. The fits by the improved GOA 
and the topGOA are shown as curves. Again the topGOA 
is superior at the moderate deformation of this example. 

We now carry out the integrals to get the angular mo- 
mentum projected energies 



E qq>tJ=0 = 2I - 1 (q,q') dcos(8) (q\HRg\q'). (6) 
Jo 
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FIG. 3: Hamiltonian matrix element as a function of rota- 
tion angle for 208 Pb at q — —500 fm 2 . The circles show the 
calculated values with the full many-particle wave function, 
and the lines are fits to the forms ho — h? sin 2 (#) (solid) and 
ho — h^il.Q — cos 6*) (dashed). The two filled circles show the 
matrix elements that were used for the fit. 



We report here the diagonal correlation energies AE qpro j 
defined by the difference in the projected energy and the 
energy of the intrinsic configuration, 



AE, 



q,pro] 



E qq j =Q 



(q\H\q). 



Later, we shall present the energy difference with respect 
to the minimum energy of the unprojected configura- 
tions, 



AE corr = mining, j =0 - whx{q\H\q). 



(7) 



This correlation energy is more useful to assess the effect 
of the projection on the calculation of the ground state 
binding energy. 

For the case presented in Figs. and 01 the energy 
gained by projection is 5.8 MeV for the accurate integra- 
tion, and essentially the same value for the GOA. In this 
case, both the topGOA and the naive GOA give values 
within 0.1 MeV of each other. In Fig. 0] we show the 
comparison of the (numerically converged) 12-point pro- 
jection with the GOA for a number of other nuclei and 
deformations. 

The correlation energy associated with the angular mo- 
mentum projection has a range of about 2 to 4 MeV, with 
the larger values occurring for well-deformed configura- 
tions. The topGOA tracks these values closely, with a 
typical error of the order of 100-200 keV. 

In the above, we only examined overlaps and matrix 
elements between configurations with the same sign of 
deformation. These quantities are of course also needed 
for matrix elements between oblate and prolate configu- 
rations. However, due to the choice of the same symme- 
try axis for oblate and prolate configurations, the overlap 
function has a rather different appearance. An example 
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FIG. 4: Comparison of explicit and GOA calculations of the 
projection, for the nuclei 120 Sn and 156 Sm and a number of 
deformations. 
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FIG. 5: Same as Fig. [2] for the off-diagonal overlap [q\R$\q') 
with q — 500 fin 2 and q = —500 fin 2 , which corresponds to 
\/h\ = 0.063. 



is shown in Fig. |SJ the overlap between the configura- 
tions q — 500 and —500 fm 2 of the nucleus 208 Pb. Here 
we see that the overlap peaks at 90°. This is not sur- 
prising. The maximum density overlap is obtained by 
aligning the long axis of the prolate ellipsoid to one of 
the longer axes of the oblate ellipsoid. This requires ro- 
tating one of the ellipsoids by 90°. We may continue to 
use the functional form of the topGOA in this case, but 
the first point should be taken where the overlap is max- 
imum. Thus, two rotations are needed, a first one by 90° 
for the maximum, which is equivalent to a permutation of 
axes, and at some intermediate angle to get the width of 
the Gaussian. The curve in Fig. [S] shows the fit obtained 
using the two points shown by filled circles. 

We determine the two parameters of the Hamiltonian 
matrix elements, eq. (JjjJ, using the same two angles. The 
fit is compared to the numerically calculated values in 
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FIG. 6: Same as Fig. |3] for the off-diagonal Hamiltonian 
matrix element {q\H Re\q') with q = 500 fm 2 and q' = —500 
fm 2 . 

Fig. H3 We can define a correlation energy for the off- 
diagonal matrix elements similar to the definition for di- 
agonal elements, 

M>r<y MV=0 <<zK /2 |-<z') ' 

For the case presented in Figs. [5] and |HJ the value cal- 
culated with the fine angular mesh is 2.50 MeV, to be 
compared with 2.49 MeV estimated by the topGOA. For 
crossover matrix elements such as this one, one should 
not expand about 9 — using a nontopological GOA. 

It is possible to use one of the other GOA's if one makes 
a rotation of one of the intrinsic configurations by 90° to 
align the longer axes. Calling the rotated wave function 
\q,n/2) = Rn/2\q), the overlap integral becomes 

I a (q,q')=2f dsm(9)( q ,TT/2\Rg\q'). (8) 
Jo 

In this representation, the improved GOA gives a value 
for the overlap which differs by less than 1% from the 
value obtained by 12-point Gauss-Legendre integration. 
Note that topGOA gives identical values with either eq. 

© or inj- 
ur MIXING DEFORMATIONS 

In the previous section, we have presented a numerical 
approximation which allows one to calculate the matrix 
element projected on J = between any two intrinsic 
configurations (q,q'). If the total number of configura- 
tions is N, one needs N(N + l)/2 elements per matrix. 
In this section, we introduce an approximation scheme 
which requires only that diagonal and nearest-neighbor 
off-diagonal matrix elements be calculated with the full 
many-particle wave functions. This reduces the number 
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of needed elements to 2N — 1, giving considerable savings 
in computational cost. 

The diagonalization of the discretized Hill- Wheeler 
equations is equivalent to the variational theory based 
on a wave function of the form 



|*) =$3ci|?i,0), 



where the Cj are variational parameters and \qi, 0) are the 
J-projected intrinsic configurations. The total energy is 
thus 



Erw = rnin ■ 



and the energy gained by configuration mixing is 



AEhw = Ehw - xmaE qq j=o- 
q 



(9) 



As mentioned earlier, the configurational basis only in- 
cludes axially symmetric deformations. For a given sign 
of the deformation (prolate or oblate), the intrinsic states 
can be ordered with respect to their quadrupole mo- 
ments. The diagonal and nearest-neighbor off-diagonal 
elements must be calculated from the full mean-field wave 
functions, but the remaining matrix elements can be es- 
timated using the GO A as follows. First consider the 
overlaps. The number and angular-momentum projected 
wave functions are renormalized to unity in the formulas 
below, i.e. we take IoiliTli) = 1- The overlaps between 
neighborin g st ates 7o(?i> <7i+i) are then used to define a 
variable x 



ti-i + A /-21og(|/ (g.. 



i,?t+l)|). 



(10) 



This plays the role of a coordinate for the GOA. Ac- 
cordingly, the overlaps between more distant states are 
estimated as 



Io(qi,qj) ~ exp(-(aii - Xj) 2 /2). 



(11) 



Let us see how well this works for the configuration set 
used for 120 Sn. The set contains 13 states, with deforma- 
tions covering the range from strongly oblate (#=-1000 
fan 2 ) to strongly prolate (g=1300 fm 2 ). A typical sep- 
aration between neighboring states is Xi+\ — x% ~ 0.8. 
Fig. [7] shows the overlaps between the second-nearest 
neighbors, plotted as an equivalent separation y\ = 

y— 21og(|/o(<Z-i, 9i+i)D- The x-axis gives the separation 
as determined from the assigned x values from eq. I)10|l. 
Ax — Xi + i — Xi-i. The accuracy of this GOA can be 
judged by the deviation of the points from the diagonal 
line. We see that eq. (jl 1|> gives satisfactory results except 
for one point. The bad point corresponds to the oblate- 
to-prolate matrix element Jo (200, —200); all other points 
correspond to transitions between states have the same 
sign of deformation. The prolate-oblate matrix element 
has a value very close to unity, I o (200, —200) = 0.996 and 
j/ 2 ' = 0.08, despite the fact that the overlaps with the 
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FIG. 7: Overlaps of angular-momentum and particle-number 
projected states from second-nearest neighbor configurations 
of 120 Sn. The ordinate and abscissa indicate the actual and 
the GOA-estimated overlaps, respectively. 



spherical state that separates them is smaller (Ax = 0.17 
and 0.19, for the prolate and oblate states, respectively). 

Unexpectedly large overlaps between prolate and 
oblate projected states were also found in the previ- 
ously mentioned study of the interacting boson model 
[25|. In fact it is easy to understand this behavior. For 
small deformations, the intrinsic configuration may be 
expanded as\q) « |0)+aQ t |0)+a 2 Q t Q t |0)/2 + . .. where 
Q is a particle-hole operator transforming as J = 2, 
M = under rotations and a is a deformation pa- 
rameter. On projection, the wave function becomes 
Pj=o\q) « |0)+a 2 /2v / 5(O t Q t )°|0) + .... This wave func- 
tion does not depend on the sign of a. Thus, the same 
state can be obtained by projecting a prolate or an oblate 
intrinsic state. 

This behavior is also found for crossover matrix ele- 
ments between more highly deformed configurations, and 
it prevents us from defining a linear metric x to estimate 
the crossover matrix elements. A possible procedure 
could be to carefully select the sequence of deformations 
so that eq. 111|) is not applied to estimate high-overlap 
crossover matrix elements. This strategy is illustrated 
in Fig. |S| The dashed line shows the path for a naive 
scheme for assigning a distance between configurations. 

The points marked A and B would be separated by 
a distance Ax given by the sum of the distances be- 
tween the origin and those two points. But from what we 
said before, the angular-momentum projected configura- 
tions A and B would actually have a very high overlap, 
larger than with the spherical configuration. Modifying 
the path by taking the short-cut AB (dotted line), the 
overlap Io(qA, <7b) will be determined exactly and the er- 
rors for more distant crossover overlaps will be reduced. 
This procedure offers an improvement, but more distance 
crossover overlaps such as AC are still underestimated. 
This suggests that the two-dimensional geometry of Fig. 
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FIG. 8: Linear metrics for GOA mixing prolate and oblate 
shapes. The naive metric is defined by the path passing 
through the spherical point, shown by the dashed line in the 
figure. An improved metric bypasses the spherical point via 
the path AB. 

[8] might be used to define a better metric for the over- 
laps: measure the distance along the lines for the pro- 
late and oblate configurations separately as done before, 
but use the triangular geometry of Fig. |H1 to determine 
crossover overlaps. We put the prolate configurations 
along the horizontal axis, and write their position as a 
vector Xi = (xi , 0) . The oblate configurations are put 
on a line making an angle a with respect to the hor- 
izontal axis, located at positions Xi = Xi (cos a, sin a). 
The distance between configurations qi, qj is taken to be 
Uij = \xi — Xj\, which is to be applied in general without 
distinguishing prolate and oblate. We will call this the 
triangular GOA, and will apply it below. In general, the 
angle IA0B of the triangle depends on the distances, be- 
coming more acute the closer A and B are to the spherical 
point. For the configuration sets used here, the spacings 
are such that the angles range from 30° to 60°. In the 
application below, we assume a fixed angle of 45°. 

We now turn to the estimation of the Hamiltonian ma- 
trix element. The first step for a GOA estimate is to ex- 
press the matrix element as a product of the overlap and 
a smoothly varying factor h{qi, qj), 

(qi\H\qj) = h(qi,qj)Io(qi,qj). 

The factor h is expanded as a power series in \xi — Xj\, 

h(qi, qj) = h {x) - h 2 (x)\xi - Xj\ 2 , (12) 

where x = (xi + Xj)/2. This equation will be used to de- 
fine h 2 between nearest-neighbor configurations qi,qi+i, 
taking the computed Hamiltonian matrix element as in- 
put. To assign nearest neighbors, we drop the spherical 
point and use the dotted path shown in Fig. |5J The def- 
inition of h 2 on these links is thus 

h 2 (ii+l) = 1 / h(qi,qi) + h(q i+1 ,qi +1 ) 

-h(qi,q i+1 )\. 

We average the h 2 on the links between two states to get 
an estimated h 2 to use in their matrix element. Then the 



matrix element of the Hamiltonian on the J— projected 
states is estimated as 

H<li,<lj) = 2 [M?»> ft) + Klh + h 2 \Si - Xj\ 2 . (13) 

We now have complete matrices both for the overlap 
and the Hamiltonian and it is diagonalized exactly the 
same way as in the full theory. This is usually done 
by first diagonalizing the matrix of overlaps Iq and con- 
structing its square root. Then the matrix 

.-1/2 o-r-1/2 

is then diagonalized, and its eigenvalues give the ener- 
gies of the GCM theory. We remind the reader that 
the method has numerical instabilities if the configura- 
tions are too close. This is circumvented by truncating 
the basis, removing states whose norm is small (typically 
< 0.01). In such situations, there is a sensitivity to el- 
ements of the norm matrix that are quite removed from 
the main diagonal. It is for this reason that the GOA 
needs to be rather accurate. As an example, the config- 
uration set used to represent 120 Sn, three of the thirteen 
eigenstates of the norm matrix had eigenvalues less that 
0.01 and were discarded. 

We now compute the correlation energy for a sample 
of nuclei using the triangular GOA and compare with 
the results of the fully calculated matrices. The resulting 
configurational correlation energy eq. © is displayed in 
Fig. OH The horizontal axis shows the value in the GOA 
using eqs. (|ll|l and JT3J, and the vertical axis shows the 
accurate result of the full calculation. 

The errors using the estimated matrix elements are all 
within the targeted accuracy of ±0.2 MeV. We conclude 
that the procedure has sufficient accuracy to be used in 
a first survey of the nuclear mass table. 

For numerical purposes in general, it is also impor- 
tant to know how close together the configurations need 




0.0 0.2 0.4 0.6 0.8 1.0 1.2 
AE RW (GOA) 

FIG. 9: Comparison of configurational correlation energies 
AEhw, eq. ©, computed with the triangular GOA to the 
numerically accurate values. 
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to be to get reliable energies. Ref. 19] indicates that 
Sx = 1.6 (overlap of 30%) is adequate and Sx = 1.1 
(overlap of 55%) is safe, taking a criterium of 0.2 MeV 
for energy convergence. This seems to be valid for the 
J— projected configurations as well. For example, for the 
nucleus 120 Sn, we can thin the configuration set taking 
out states whose overlaps with neighbors exceeds 55%. 
Here the thinned space has only 4 configurations, com- 
pared with 13 in the original set. The configurational 
correlation energy is nearly unchanged by the thinning 
and is within the desired error bound. Clearly, more 
study is needed to find an optimum spacing for a global 
prescription. 



IV. EXAMPLES OF CORRELATION ENERGIES 

In this section we survey the trends of the correla- 
tion energy taking the results for a number of nuclei. In 
this sampling we have chosen nuclei to represent differ- 
ent kinds of structure in the mean field approximation, 
including closed-shell magic ( 208 Pb), strongly deformed 
( 156 Sm), and semi-magic with strong pairing ( 120 Sn). In 
addition, we consider two isotope chains, krypton and 
lead, which exhibit coexisting oblate and prolate defor- 
mations. The results are presented in Table [I] The 
first three columns of energies show the numerically con- 
verged correlation energies obtained with the N, Z, and 
J projections and Hill- Wheeler representation as in refs. 

mum 

The first column gives AE corr , the energy gain over 
that of the lowest intrinsic configuration associated with 
angular momentum projection, defined in eq. 0. We 
should note that we only use the given configurations to 
determine the energy minimum for the projected states, 
min 9 E qq j=o, which could give some error in assigning 
the energy gain. There will be a compensating error in 
the value of the configurational correlation energy, with 
the total independent of min g E qq ^j = Q. In the table, the 
nuclei for which the deformation q is significantly differ- 
ent for the intrinsic and the projected configurations are 
indicated with asterisks. An example where the deforma- 
tion is virtually identical before and after projection is 
156 Sm. Its correlation energy is 3.0 MeV, which is rather 
typical for deformed nuclei. What is surprising to us is 
that the other nuclei have rotational correlation energies 
of the same order of magnitude. The semi-magic "spheri- 
cal" nucleus 120 Sn is a case in point. At the mean-field or 
intrinsic level, the ground state is indeed spherical, but 
the energy after angular momentum projection is lower 
for the configuration with q = 400 fm 2 . Qualitatively 
this behavior is to be expected, because the variation- 
after-projection is a way to introduce correlations from 
the attractive parts of the interaction that only have per- 
turbative character. The surprising finding is the mag- 
nitude, which is only 10% lower in 120 Sn than the value 
for a fully deformed nucleus. In fact all the nuclei except 
the doubly-magic 208 Pb have correlation energies in the 



TABLE I: Correlation energies in the GCM associated with 
quadrupole deformations. Energies are given in MeV. 



z 


A 


/ \ Hi 


/\ Hi £j~\Xr 


total GOA 


36 


74 


2.9 


0.8 


3.7 


3.7 


36 


78 


2.4* 


0.9 


3.3 


3.2 


36 


82 


3.0 


0.8 


3.8 


3.5 


36 


86 


2.8 


0.4 


3.3 


3.0 


50 


120 


2.8* 


0.8 


3.6 


3.3 


62 


156 


3.0 


0.8 


3.8 


3.8 


82 


186 


2.7 


0.9 


3.6 


3.3 


82 


190 


2.7 


1.0 


3.8 


3.7 


82 


194 


2.8 


1.1 


3.9 


3.7 


82 


208 


1.5* 


0.2 


1.7 


1.6 



same range. 

The next column of the table shows additional energy 
gain obtained by solving the Hill- Wheeler equation to 
mix configurations, AEhw i n eq. ©■ Except for the 
doubly-magic 208 Pb, these numbers are remarkably con- 
stant over the nuclei considered, ranging in value from 
0.8 to 1.1 MeV. The nucleus 208 Pb behaves in a simi- 
lar way as the model discussed in ref. (2^. That model 
permitted a significant energy by projection, but the op- 
timal projected state had a very high overlap with the 
true eigenstate, and there was no further energy gain by 
mixing configurations. 

The third column shows the total correlation energy, 
AE corr + AEhw- The range of variation is 0.6 MeV, 
except for 208 Pb. 

The last column in Table [I] shows the total correlation 
energy computed with the approximations in the treat- 
ment of the Hill- Wheeler equation described in the last 
section. The results are very close to that of the full treat- 
ment, with the r.m.s. error having a value of 0.15 MeV. 
We now come to the conclusions that we draw from these 
results. 



V. CONCLUSIONS 

Although the primary purpose of this study was nu- 
merical and we only examined a small sample of nuclei, 
our findings suggest the possible physics that may emerge 
from a systematic, global study of the quadrupolar corre- 
lation energy. One might naively expect that the energy 
of projection would be large for deformed nuclei but not 
for others. In fact we found that this energy is rather 
large except for the doubly magic nucleus that we con- 
sidered. The additional energy associated with vibrations 
is small, less than one MeV, and not very different from 
one nucleus to another. For the sum we see fluctuations 
of the order a few hundred keV. 

In the global theories of nuclear masses, one treats the 
mean-field effects in detail but up to now ignoring or ap- 
proximating in a phenomenological way fluctuating parts 
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of the correlation energy. The relative constancy of our 
calculated correlation energy suggests that these fluctua- 
tions are indeed small, permitting a mean-field treatment 
to be quite successful. On the other hand, the fluctuation 
effects we have found are large enough to be potentially 
useful. In the present mass theories there are r.m.s. de- 
viations of the order of 700 keV. For the nine nuclei in 
Table[U the r.m.s. fluctuation in the total correlation en- 
ergy is 650 keV. This gives grounds for the hope that 
including the quadrupolar correlation energy in calculat- 
ing nuclear masses might significantly improve the agree- 
ment with the experimental values and produce a theory 
with better predictive power. 

On the numerical side, the GOA for angular momen- 
tum projection reduces the computational cost by almost 
an order of magnitude. A similar saving is realized by 
the triangular GOA for the Hill- Wheeler matrices, which 
only requires computing diagonal and next-to-diagonal 
elements. This brings the total cost down the same level 
as the HF+BCS computations of the GCM configura- 
tions, and makes feasible a systematic determination of 
correlation energies for the nuclear mass table. 

We thus have the intention of continuing this project 
to compute the quadrupolar correlation energies of all 
even-even nuclei. For a global calculation we would need 
to use an energy functional that is applicable across the 
chart of nuclides. We have been using the SLy6 Skyrme 
parameterization, but so far the pairing parameters have 
been different for light and heavy nuclei. Ultimately, a 
new fit of the mean field energy functional will be re- 
quired, but for a first survey of the correlation energy it 
should be sufficient to take the same functional, but with 
a fixed pairing interaction. 
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Appendix: Summary of computational costs 

Our aim is to reduce the computing time needed to cal- 
culate correlation energies at a level comparable to that 
of mean-field calculations. Let us estimate the cost of 
both and to determine how they scale with the number 
of nucleons. The parameters which govern the comput- 
ing time are the number of active orbitals N act and the 
dimension of the vector representing an orbital N r . The 



mesh size is always fixed to 0.8 fm, a value sufficient to 
have a satisfactory accuracy on energies. The number of 
mesh points must be large enough to guarantee that the 
wave functions can be set equal to zero at the edge of 
the box so that derivatives and the Coulomb boundary 
conditions can be calculated accurately. The mean-field 
equations are solved by the imaginary time method and 
one can limit the number of orbitals that have to be com- 
puted to the orbitals close to the Fermi level. In practical 
applications, to be on the safe side, N act and N r are cho- 
sen larger than required. 

In Hartree-Fock theory N act is equal to N occ , the num- 
ber of occupied orbitals. Due to time-reversal symmetry, 
the orbitals are two-fold degenerate and the number of 
orbitals is half the number of nucleons. When pairing 
correlations are taken into account, either at the BCS or 
at the Bogoliubov level, N act is larger. In practical ap- 
plications, to be sure that enough HF wave functions are 
computed, N occ is taken around 2A/Z. 

Calculations are performed in a 3-dimcnsional coordi- 
nate grid representation. To calculate the nuclear ground 
state and to be able to easily rotate the orbital wave func- 
tions on the mesh, we consider only cubic meshes. In 
practice, the actual dimension of the grid can be reduced 
by symmetry considerations. The reflection symmetry of 
the mean-field Hamiltonian allows one to consider only 
the points in an octant. We will thus define N x to be 
the number of points along the positive x-axis. How- 
ever, there is a price to pay for using this symmetry: the 
wave functions must be complex. Also, the spin degree 
of freedom doubles the size of the vector. In the end the 
dimension of a real array representing an orbital is 

N r = AN 3 x . 

The number of points scale roughly as 3A 1 / 3 , although 
the number used in practical applications is larger for 
light nuclei. This dependence on A ensures that the wave 
functions can be safely put to zero at the edge of the box. 

The mean-field calculation is done iteratively, deter- 
mining in each cycle the total energy of the nucleus, the 
HF Hamiltonian and its action on the individual wave 
functions. These computations need to be repeated a 
number of time equal to the total number of iterations 
required for a full convergence, typically 200 to 400 times. 
It has also to be done for each quadrupole moment which 
will enter the GCM calculation. The GCM requires cal- 
culating the overlaps and the Hamiltonian matrix ele- 
ments between wave functions with different orientations 
and deformations. Let us first evaluate the computa- 
tional effort that is required for a single set of Euler an- 
gles and a single deformation. This effort is similar to 
that required to calculate the energy in the mean-field 
calculation. 

For the calculation of a matrix element of the GCM 
Hamiltonian kernel or the total HF energy the dominant 
terms are those with derivatives of the wave functions. 
This has a computational cost given by 

6N act N r N x , (14) 
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following the method of ref . [2lJ ■ A prefactor 3 accounts 
for the three Cartesian directions along which derivatives 
have to be calculated. The prefactor 2 counts separately 
the multiplication of the initial vector and the addition to 
the final vector. Several terms in the Hamiltonian matrix 
element involve derivations, e.g. those associated with 
the spin-orbit field, and the total computational effort is 
an order of magnitude larger than given by the expression 
Eq. (|14fl . More terms appear in the GCM calculation due 
to the presence of vector densities which are equal to zero 
for diagonal matrix elements. 

The cost related to particle-number projection is not 
easy to evaluate, since this projection is intimately mixed 
with the other part of the calculation. We use the method 
of ref. Q which requires one to multiply the v factor of 
the Bogoliubov transformation by a phase. The calcu- 
lation of the derivatives of the wave functions is not af- 
fected by this projection. However, all densities become 
complex. 

The mean-field calculation requires one to calculate 
also the action of the HF Hamiltonian on all mean-field 
wave functions. Once more, the most time-consuming 
part concerns terms including derivatives. They are 
grouped in terms with first and second order derivatives 
acting on the four components of the wave functions, with 
24 terms in total. 

In the GCM case, another time-consuming task is ro- 
tating the many-particle wave function, an operation that 
is needed for angular momentum projection. This is done 
by interpolating the single-particle wave functions in the 
x — z plane to the points corresponding to some rotation 
9 about the y-axis, 

Rg : (x, z) — (x — > x cos 8+z sin#, z — » ~x sin 8+z cos 9). 

The needed accuracy is achieved by using all the points 
in a half plane to interpolate in that plane. The compu- 
tational cost here is 

AN act N r N 2 x 

where the factor 4 includes a factor 2 for expanding from 
a quarter-plane slice of the octant to the half plane. 



To take overlaps, one needs to calculate a determinant. 
The determinant is of order N act and setting up its ma- 
trix requires 

2/V 2 N 



operations. Since N r is usually less than N act , this is 
usually the most time consuming task. However, the si- 
multaneous projection on particle number and angular 
momentum requires also additional calculations to follow 
the phase of the overlap kernel. Test calculations indi- 
cate that the particle number projection approximately 
doubles the computing of the determination of the GCM 
kernels. 

For our calculations, N x is of the order of 14 (Kr) to 
16 (Pb). Thus the vector size is about N r s» 4 * 15 3 = 
13, 500. The number of orbitals for Pb is around N act — 
140. With these numbers, the interpolation of the wave 
functions has about the same cost as setting up the de- 
terminants. 

Putting all these numbers together with realistic val- 
ues, one can estimate that the determination of a single 
point for the GCM+projection calculation takes a factor 
10 longer than a single mean-field iteration. 

The computing time is also determined by the number 
of times these basic calculations have to be done. It is the 
number of iterations (typically 300) times the number of 
quadrupole moments that have to be considered (typi- 
cally 15) for the mean-field part. For the GCM part, it 
is the number of active Euler angles (typically 12) times 
the number of matrix elements (typically 15 2 ). With all 
these numbers, the GCM calculation is an order of mag- 
nitude longer than the mean-field calculation. Reducing 
the number of Euler angles to 1 or 2 and the number of 
matrix elements to be determined exactly to 2 times the 
number of discretized quadrupole moments reduce this 
part of the calculation to a time similar to the mean-field 
calculation. 
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